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Abstract 

This paper provides an in-depth empirical analysis of several evolutionary algorithms on the 
one-dimensional spin glass model with power-law interactions. The considered spin glass model 
provides a mechanism for tuning the effective range of interactions, what makes the problem 
interesting as an algorithm benchmark. As algorithms, the paper considers the genetic algorithm 
(GA) with twopoint and uniform crossover, and the hierarchical Bayesian optimization algorithm 
(hBOA). hBOA is shown to outperform both variants of GA, whereas GA with uniform crossover 
is shown to perform worst. The differences between the compared algorithms become more 
significant as the problem size grows and as the range of interactions decreases. Unlike for GA 
with uniform crossover, for hBOA and GA with twopoint crossover, instances with short-range 
interactions are shown to be easier. The paper also points out interesting avenues for future 
research. 

Keywords: Spin glass, power-law interactions, hierarchical BOA, genetic algorithm, estimation of 
distribution algorithms, evolutionary computation. 

1 Introduction 

Spin glasses are prototypical models for disordered systems, which provide a rich source of difficult 
computational problems. Prom the optimization perspective, the problem of finding ground states 
of various spin glass models represents a popular class of challenging optimization problems useful 
for analysis, design, and enhancement of optimization algorithms. There are a variety of spin glass 
models available, from short-range models arranged on a 2D lattice to the infinitely dimensional 
Sherrington-Kirkpatrick (SK) model [23] . One of the interesting models is the one-dimensional spin 
glass model with power-law interactions [251 [9] where the range of interactions can be controlled via 
a parameter. This provides the user with enough flexibility to cover a broad spectrum of spin glass 
models from short-range models, such as the one-dimensional nearest-neighbor spin glass, to long- 
range models, such as the infinitely dimensional SK spin glass. The one-dimensional model with 
power-law interactions has been recently used to elucidate some spin glass properties and suggested 
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as an interesting algorithm benchmark [22\ \2~T\ I20| . Although spin glasses have been extensively 
used as test problems for evolutionary algorithms, most studies focused on short-range and finite 
dimensional systems such as the 2D or 3D spin glass [281 E21 US US [38l E HI]- Furthermore, the 
one-dimensional spin glass with power-law interactions has not yet been studied in the context of 
evolutionary algorithms. 

This paper presents an in-depth empirical analysis of several evolutionary algorithms on the one- 
dimensional spin glass model with power-law interactions. In the comparison, the paper considers 
the genetic algorithm with two crossover operators and the hierarchical Bayesian optimization 
algorithm. The effects of the parameter that controls the range of interactions are studied and for 
each setting of this parameter a range of problem sizes are analyzed. The results are analyzed and 
related to known properties of the evolutionary algorithms under consideration. Important topics 
for future work in this area are also outlined. 

The paper is organized as follows. Section [2] describes the one-dimensional spin glass model with 
power-law interactions. Section [3] outlines evolutionary algorithms included in this comparison. 
Section 2] presents and discusses experimental results. Section [5] focuses on interesting topics for 
future work in this area. Finally, section [6] summarizes and concludes the paper. 

2 One-Dimensional Spin Glass with Power-Law Interactions 

The Sherrington-Kirkpatrick spin glass [23] is described by a set of spins {s^} and a set of couplings 
{Ji,j} between all pairs of spins. For the classical Ising model, each spin Sj can be in one of 
two states: s, = +1 or Sj = —1. This simplification corresponds to highly anisotropic systems; 
nevertheless, the two-state Ising model comprises all basic effects also found in models with more 
degrees of freedom. 

For a set of coupling constants {Jj and a configuration of spins C = the energy can be 
computed as 



There are two typical tasks in statistical physics: (1) integrate a known function over all pos- 
sible configurations of spins for given coupling constants assuming the Boltzmann distribution of 
spin configurations; (2) find one ore more ground states (spin configurations associated with the 
minimum possible energy). The problem of finding ground states is NP-complete even when the 
interactions are limited only to neighbors in a 3D lattice [2j; the SK spin glass is thus certainly 
NP-complete unless we severely restrict couplings. 

Typically, one would study a number of randomly generated spin glass instances where couplings 
for each instance would be generated according to the Gaussian distribution with zero mean and 
unit variance. In the one-dimensional spin glass model with power-law interactions [20], spins are 
arranged equidistantly on a circle with circumference n. The spins are numbered in a counterclock- 
wise fashion. Every spin interacts with every other spin, just like in the standard SK spin glass; 
however, interactions between spins located further from each other are weaker, and the effects of 
distance can be controlled with a parameter, providing a mechanism for tuning the effective range 
of interactions. More specifically, the couplings are generated according to 
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Figure 1: One-dimensional spin glass of size n = 10 arranged on a ring. 



where e^j are generated according to normal distribution with zero mean and unit variance, c(cr) 
is a normalization constant, a > is the user-specified parameter to control the effective range of 
interactions, and rjj = nsin(7r|i — j\/n)/ir denotes the geometric distance between S{ and Sj (see 
figure [1]). The magnitude of spin-spin couplings decreases with their distance. Furthermore, as 
discussed shortly, the effects of distance on the magnitude of couplings increase with a. 

By varying a one can tune the model from the infinite-range to the short-range universality 
class: For < a < 1/2 the model is in the infinite-range universality in the sense that X^I^W 
diverges, and for a = it corresponds to the SK model. For 1/2 < a < 2/3 the model describes 
a mean-field long-range spin glass, corresponding — in the analogy with short-range systems — to a 
short-range model in dimension above the upper critical dimension d > d u = 6. For 2/3 < a < 1 
the model has non-mean-field critical behavior with a finite transition temperature T c , and for 
a > 1, the transition temperature is zero and the behavior of the model is short ranged. 

To find guaranteed ground states, a branch-and-bound algorithm adopted from refs. |17[ 123] 
was used. This allows us to find guaranteed ground states of instances of up to about n = 100 spins 
for larger values of a and up to about n = 80 spins for smaller values of a. For larger systems, we 
used the population-doubling approach proposed in ref. [32]. 

3 Compared Algorithms 

The genetic algorithm (GA) [THl [UJ evolves a population of candidate solutions typically repre- 
sented by binary strings of fixed length with the first population generated at random according 
to the uniform distribution over all binary strings. Each iteration starts by selecting promising 
solutions from the current population; we use binary tournament selection without replacement. 
New solutions are created by applying variation operators to the population of selected solutions. 
Specifically, crossover is used to exchange bits and pieces between pairs of candidate solutions and 
mutation is used to perturb the resulting solutions. Here we use uniform or twopoint crossover, and 
bit-flip mutation |11| . To maintain useful diversity in the population, the new candidate solutions 
are incorporated into the original population using restricted tournament selection (RTS) [T3J. The 
run is terminated when termination criteria are met. In this paper, each run is terminated ei- 
ther when the global optimum has been found or when a maximum number of iterations has been 
reached. 

The hierarchical Bayesian optimization algorithm (hBOA) [371 [39] [36] is an estimation of dis- 
tribution algorithm (EDA) [U EQ] [26] SO] [27] [43]. EDAs— also called probabilistic model-building 
genetic algorithms (PMBGAs) [30] and iterated density estimation algorithms (IDEAs) [6] — differ 
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from GAs by replacing standard variation operators of GAs such as crossover and mutation by 
building a probabilistic model of promising solutions and sampling the built model to generate 
new candidate solutions. The only difference between GA and hBOA variants used in this study 
is that instead of using crossover and mutation to create new candidate solutions, hBOA learns a 
Bayesian network with local structures [7J [10] as a model of the selected solutions and generates 
new candidate solutions from the distribution encoded by this model. For more details on hBOA, 
see refs. (37J 136] . 

The deterministic hill climber (DHC) is incorporated into both GA and hBOA to improve their 
performance similarly as in previous studies on using GA and hBOA for solving spin glasses [38U36] . 
DHC takes a candidate solution represented by an n-bit binary string on input. Then, it performs 
one-bit changes on the solution that lead to the maximum improvement of solution quality. DHC is 
terminated when no single-bit flip improves solution quality and the solution is thus locally optimal. 
Here, DHC is used to improve every solution in the population before the evaluation is performed. 

To represent a configurations for n spins, both hBOA and GA use n-bit binary strings with 
each bit corresponding to one of the spins; state +1 is represented by a 1, state —1 is represented 
by a 0. 

4 Experiments 

4.1 Problem Instances 

Instances of sizes n = 20 to n = 150 were considered, and for each problem size, the coefficient a 
was set to each of the values from {0.00, 0.55, 0.75, 0.83, 1.00, 1.50, 2.00} (the set of values of a was 
adopted from ref. [5]). For each combination of n and a, 10,000 problem instances were generated 
and tested. This resulted in the total of 610,000 unique problem instances. Instances of up to 
n = 100 were solved using the branch and bound algorithm; other instances were solved with the 
population-doubling method and hBOA as described in ref. [32] . 

4.2 Experimental Methodology 

For each problem instance and each algorithm, an adequate population size was approximated using 
bisection [46 \ [36] to ensure that the optimum is found in 10 out of 10 independent runs. Each run 
was terminated when the global optimum was found (success) or when the maximum number of 
generations n was reached (failure). The most important statistic relating to the overall complexity 
of each algorithm is the number of DHC flips, since this statistic combines all important statistics 
and can be consistently compared regardless of the used algorithm. Due to the variety of hardware 
used to perform the simulations, it was not possible to compare the actual running times. 

4.3 Results 

The growth of the number of evaluations and the number of DHC flips until optimum for hBOA, GA 
with twopoint crossover, and GA with uniform crossover is shown in figures [2] and [3l As expected, 
both the number of evaluations as well as the number of flips grow with problem size, and a 
similar observation can be made for the adequate population size (see figure [3J and the number 
of generations (results omitted). This result is not surprising because the more decision variables 
in the problem, the more difficult the problem becomes. What is somewhat more interesting is 
how the rate of growth of these statistics changes with a, which controls the range of interactions, 
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Figure 2: Growth of the number of evaluations with problem size. 
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Figure 3: Growth of the number of flips with problem size. 



and how the effects of a change depending on the algorithm under consideration; this is the topic 
discussed in the following few paragraphs. 

Based on the definition of the ID spin glass with power-law interactions, as the value of a grows, 
the range of the most significant interactions is reduced. With reduction of the range of interactions, 
the problem should become easier both for selectorecombinative GAs capable of linkage learning, 
such as hBOA, as well as for selectorecombinative GAs which rarely break interactions between 
closely located bits, such as GA with twopoint crossover. This is clearly demonstrated by the 
results for these two algorithms presented in figures [2] and [3l Although for many problem sizes, the 
absolute number of evaluations and the number of flips are in fact smaller for larger values of a, 
the growth of the two statistics with problem size slows down substantially as a is increased. In 
fact, for the smallest values of a, the number of evaluations and the number of flips both appear to 
grow faster than polynomially, whereas for the largest values of a, the growth appears to be upper 
bounded by a polynomial of low order. On the other hand, for the uniform crossover, which does 
not respect interactions between consequent bits, the growth appears to be faster than polynomial 
for both the number of evaluations and the number of flips. 

To better visualize the effects of a on the performance of compared algorithms, figures [5] and [6] 
show the ratio of the number of evaluations and the number of flips for a < 2 and a = 2. The 
larger the ratio, the better the results for a = 2 compared to smaller values of a. For both hBOA 
and GA with twopoint crossover, the ratios for a < 1 appear to grow relatively fast; in fact, the 
growth appears to be faster than polynomial; on the other hand, for a = 1.5, the growth can still 
be observed but it is very slow in both cases. For GA with uniform crossover, the opposite behavior 
can be observed; this indicates that uniform crossover does not benefit from the decrease in the 
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Figure 4: Growth of the population size with problem size. 
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Figure 5: Comparison of the number of evaluations for a = 2.00 with that for a < 2.00. 



range of interactions, at least within the range of values considered in this work. 

To compare performance of two algorithms on any specific problem instance, one may compute 
the ratio of the number of evaluations required by the two algorithms; analogical comparison can 
be performed for the number of flips, the population size, or the number of generations. These 
ratios can then be averaged over all instances for any considered combination of values of n and a. 
The comparison for the number of evaluations via such ratios is shown in figure [7] and analogical 
results for the number of flips and the required population sizes are shown in figures [8] and 
respectively. Pairwise comparisons of hBOA and GA show that with respect to both the number of 
evaluations as well as the number of flips, the best performance is achieved by hBOA, the second 
best performance is achieved by GA with twopoint crossover, and the worst performance is achieved 
by GA with uniform crossover. The differences between all algorithms increase with problem size 
and the differences for larger values of a are more significant than those for smaller values of a; in 
other words, the shorter the range of interactions, the more significant the differences. 

One of the interesting questions from the perspective of problem difficulty is whether the prob- 
lem instances that are difficult for one algorithm are also difficult for other algorithms. One of the 
approaches to visualizing this is to look at the correlation between the running times between pairs 
of algorithms over problem instances for specific values of n and a. Figure PTOl compares the running 
times measured by the number of flips for the largest problem of n = 150 spins and a G {0.55, 2.00}. 
The results indicate, both numerically and visually, the the running times for a = 0.55 are more 
strongly correlated than the running times for a = 2.00. A somewhat more detailed view on the 
correlation between the numbers of flips until optimum for the different algorithms is shown in 
figure [Til which confirms that the results are more strongly correlated for smaller values of a; ad- 
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(a) hBOA (b) GA (twopoint) (c) GA (uniform) 

Figure 6: Comparison of the number of flips for a = 2.00 with that for a < 2.00. 
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Figure 7: Ratio for the number of evaluations for pairs of compared algorithms. 

ditionally, these results indicate that the correlations are almost consistently stronger for smaller 
problem sizes. 

5 Future Work 

The numerous statistics provided by the experiments presented in this paper can be analyzed 
in more detail to provide a deeper insight into the strengths and weaknesses of the compared 
evolutionary algorithms, and to provide inputs for further improvement of these algorithms on the 
one-dimensional spin glass with power-law interactions and on other difficult classes of additively 
decomposable problems with complex fitness landscape. 

One of the important outputs of this paper was the large number of problem instances of the 
one-dimensional spin glass with power-law interactions. These instances can be used for additional 
empirical studies in evolutionary computation and beyond. Most importantly, one should consider 
techniques that were shown to perform well on various classes of spin glass models and compare 
the performance of these techniques to that of the evolutionary algorithms studied in this paper. 
Three algorithms are of special interest: extremal optimization [3], hysteretic optimization [34^135]. 
and GA with triadic crossover [53]. Hybrids of evolutionary algorithms, hysteretic optimization, 
and extremal optimization hold a great promise for solving large instances of both the general SK 
spin glass as well as the one-dimensional SK spin glass with power-law interactions. The key is to 
identify the strengths and weaknesses of the different algorithms and to combine their strengths 
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Figure 9: Ratio for the required population sizes for pairs of compared algorithms. 



effectively. 

Another interesting topic for future research is related to problem difficulty of the one- 
dimensional spin glass with power-law interactions and the general SK spin glass. Is there a way 
to quantify features that make one problem instance more difficult than the other one besides the 
direct influence of n and er? Can scalability theory for evolutionary algorithms \13\ [T2| \15\ Wl\ and 
estimation of distribution algorithms [3lJ EHJ HU [36l EH EH US] be used to gain better understand- 
ing of problem difficulty of the considered spin glass model? Can problem difficulty be explained 
using existing techniques for landscape analysis in evolutionary computation and optimization the- 
ory [331 13 113 E] or do new tools have to be designed for this purpose? 



6 Summary and Conclusions 

The problem of finding ground states of the one-dimensional spin glass model with power-law inter- 
actions provides an interesting class of optimization benchmarks with tunable range of interactions. 
This paper presented an in-depth empirical analysis of the genetic algorithm (GA) with twopoint 
and uniform crossover and the hierarchical Bayesian optimization algorithm (hBOA) on this class 
of problems. A large number of instances of different size and range of interactions were gener- 
ated. The branch and bound algorithm was then used to find guaranteed ground states of most 
instances; for largest instances, which were unsolvable in practical time by the branch and bound, 
the population-doubling approach based on hBOA was used to find reliable (but not guaranteed) 
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optima. Overall, 610,000 unique problem instances were used. All algorithms were then applied to 
the generated instances, and the results were analyzed and discussed. 

The best performance was achieved by hBOA, the second best performance was achieved by GA 
with twopoint crossover and the worst performance was achieved by GA with uniform crossover. 
The differences between all algorithms became more substantial as the range of interactions was 
reduced and as the problem size was increased. 

The results clearly indicated that when the recombination operator is able either to learn linkage 
automatically (hBOA) or to preserve linkage between consequent bits (GA with twopoint crossover), 
instances with short-range interactions are clearly easier than those with long-range interactions. 
Nonetheless, for recombination that processes each bit independently regardless of its position in 
the solution strings (GA with uniform crossover), problem instances with long-range interactions 
are in fact easier than those with short-range interactions. 

For instances with short-range interactions, both hBOA and GA with twopoint crossover were 
empirically shown to require low-order polynomial time whether we measure time complexity 
through the number of evaluations or the number of flips of the local searcher. Nonetheless, for 
long-range interactions, scalability of both algorithms suffered and the time required to solve these 
problems appeared to grow faster than polynomially. For GA with uniform crossover, the growth 
of the time complexity regardless of the range of interactions appeared to be faster than polynomial 
for the entire range of problem instances under consideration. 

This paper presented the first empirical analysis of evolutionary algorithm on the problem of 
finding ground states of the one-dimensional spin glass with power-law interactions. We believe that 
this class of problems provides a rich source of problem instances for future research in evolutionary 
computation, both in the area of empirical testing as well as in the design of theoretical models of 
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(a) hBOA vs. GA (twopoint) (b) hBOA vs. GA (uniform) (c) GA (uniform) vs. GA (twopoint) 



Figure 11: The correlation between the number of flips for different algorithms for n = 50 and 
n = 150. 

evolutionary algorithms on additively decomposable problems with a complex fitness landscape. 
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